#################################################################################
#                         GOV 2001: REPLICATION PROJECT                         #
#  Article: Lipsmeyer, C.S. & Zhu, L. (2011)  "Immigration, Globalization, and  #
# Unemployment Benefits in Developed EU States", AJPS Vol. 55(3), Pp. 647-664   #
#                     Authors: Diana Draghici & Jerome Hughes                   #
#                       This version: Thu 22 Mar 2012                           #     
#################################################################################

#################################################################################
#       Preliminaries           
#################################################################################


# Set working directory 
setwd("//Users/DianaDraghici//Documents//HARVARD UNIVERSITY//Courses//Spring Semester 2012//GOV 2001//Replication Project//Lipsmeyer&Zhu_AJPS_2011")

# Libraries
library(foreign)
library(tseries)
library(plm)
library(AER)
library(zoo)
library(TSA)
library(pcse)
library(xtable)
library(lattice)
library(Zelig)
library(lmtest)
library(nlme)
library(mvtnorm)
library(graphics)


# Load dataset 
data <- read.dta("AJPS Immigration Data.dta")
attach(data)

### Expand the dataset to create a full set of state-year observations 
# (equivalently to -tsfill, full- in Stata);
# this is b/c there are unbalanced panels and time series of unequal length in the data,
# which makes some computations and the use of some packages impossible
tscs <- cbind(rep(1:15, each=37),  rep(1971:2007,15))
colnames(tscs) <- c("state","year")
fulldata <- merge(tscs, data, by= c("state","year"), all=TRUE)

# Convert to TSCS format & attach
tsdata <- pdata.frame(fulldata, index=c("state", "year"), drop.index=FALSE, row.names=TRUE)
attach(tsdata)



#### Data Management
### Computing lagged and first-differenced variables

  # Function to create lags
lag.fn <- function(x) {lag(x, 0:1)[,2]}
  # Apply function
tsdata$L.i_entitlement <- lag.fn(tsdata$i_entitlement)
tsdata$L.i_socialcorp3 <-  lag.fn(tsdata$i_socialcorp3)
tsdata$L.leftgov <-  lag.fn(tsdata$leftgov)

  # Function to create first-differences
diff.fn <- function(x) {diff(x, 0:1)[,2]}
  # Apply function
tsdata$D.netimmrate <- diff.fn(tsdata$netimmrate)
tsdata$D.laborcost_new <- diff.fn(tsdata$laborcost_new)
tsdata$D.unemp <- diff.fn(tsdata$unemp)
tsdata$D.trade <- diff.fn(tsdata$trade)
tsdata$D.gdp <- diff.fn(tsdata$gdp)
tsdata$D.fdi <- diff.fn(tsdata$fdi)


###############################
# Figure 1 on p.              #
##############################

# NOTE: Because the core dataset was uncoded. We did not append the legend with the country code. 
data2 <- read.dta("AJPS Immigration Data.dta")
#Generate plot 
pdf("AJPS2011_Figure1.pdf", 10, 6) 

par(mfrow=c(3,1))

plot(0, type="n", xlim=range(1970,2010), ylim=range(.6,1),
     xlab ="Year", ylab="EU Labor Market Integration Index")

lines(data2$year[state==1], data2$laborcost_new[state==1], col="black", type="l", lwd=1, lty=1)
lines(data2$year[state==2], data2$laborcost_new[state==2], col="black", type="l", lwd=1, lty=2)
lines(data2$year[state==3], data2$laborcost_new[state==3], col="black", type="l", lwd=1, lty=3)
lines(data2$year[state==4], data2$laborcost_new[state==4], col="black", type="l", lwd=1, lty=4)
lines(data2$year[state==5], data2$laborcost_new[state==5], col="black", type="l", lwd=1, lty=5)

text (1975, .18, "European Council Conclusion", cex =.6)
abline (v=1975)
text (1985, .18, "Single European Act", cex=.6)
abline (v=1985)
text (1992, .18, "Maastricht Treaty", cex=.6)
abline (v= 1992)
text (1998, .18, "Amsterdam \n Treaty", cex=.6)
abline (v=1998)
text (2002, .18, "Nice Treat", cex=.6)
abline (v=2002) 
text (2007, .18, "Lisbon Treaty", cex=.6)
abline (v=2007)

plot(0, type="n", xlim=range(1970,2010), ylim=range(.6,1),
     xlab ="Year", ylab="EU Labor Market Integration Index")

lines(data2$year[state==6], data2$laborcost_new[state==6], col="black", type="l", lwd=1, lty=1)
lines(data2$year[state==7], data2$laborcost_new[state==7], col="black", type="l", lwd=1, lty=2)
lines(data2$year[state==8], data2$laborcost_new[state==8], col="black", type="l", lwd=1, lty=3)
lines(data2$year[state==9], data2$laborcost_new[state==9], col="black", type="l", lwd=1, lty=4)
lines(data2$year[state==10], data2$laborcost_new[state==10], col="black", type="l", lwd=1, lty=5)

text (1975, .18, "European Council Conclusion", cex =.6)
abline (v=1975)
text (1985, .18, "Single European Act", cex=.6)
abline (v=1985)
text (1992, .18, "Maastricht Treaty", cex=.6)
abline (v= 1992)
text (1998, .18, "Amsterdam \n Treaty", cex=.6)
abline (v=1998)
text (2002, .18, "Nice Treat", cex=.6)
abline (v=2002) 
text (2007, .18, "Lisbon Treaty", cex=.6)
abline (v=2007)

plot(0, type="n", xlim=range(1970,2010), ylim=range(.6,1),
     xlab ="Year", ylab="EU Labor Market Integration Index")

lines(data2$year[state==11], data2$laborcost_new[state==11], col="black", type="l", lwd=1, lty=1)
lines(data2$year[state==12], data2$laborcost_new[state==12], col="black", type="l", lwd=1, lty=2)
lines(data2$year[state==13], data2$laborcost_new[state==13], col="black", type="l", lwd=1, lty=3)
lines(data2$year[state==14], data2$laborcost_new[state==14], col="black", type="l", lwd=1, lty=4)
lines(data2$year[state==15], data2$laborcost_new[state==15], col="black", type="l", lwd=1, lty=5)

text (1975, .18, "European Council Conclusion", cex =.6)
abline (v=1975)
text (1985, .18, "Single European Act", cex=.6)
abline (v=1985)
text (1992, .18, "Maastricht Treaty", cex=.6)
abline (v= 1992)
text (1998, .18, "Amsterdam \n Treaty", cex=.6)
abline (v=1998)
text (2002, .18, "Nice Treat", cex=.6)
abline (v=2002) 
text (2007, .18, "Lisbon Treaty", cex=.6)
abline (v=2007)


dev.off()

####### End of Figure 1 ############

############

### Simple OLS regression
Mod1 <- lm(i_entitlement ~ L.i_entitlement + D.netimmrate + D.laborcost_new +
  Dimmrate_Dlabcost + Dimmrate_Lsocial + Dimmrate_Lleft +
  L.i_socialcorp3 + L.leftgov + D.unemp + D.trade + D.gdp  + D.fdi +
  as.factor(state), data=tsdata)
summary(Mod1)
OLScoefs <- as.matrix(Mod1$coef)
xtable(OLScoefs, digits=4)

#######################################
# Implement Prais Winsten transformation (using John Fox's code)

prais.winsten.lm <- function(mod){
     X <- model.matrix(mod)
     y <- model.response(model.frame(mod))
     e <- residuals(mod)
     n <- length(e)
     names <- colnames(X)
     rho <- sum(e[1:(n-1)]*e[2:n])/sum(e^2)
     y <- c(y[1] * (1 - rho^2)^0.5, y[2:n] - rho * y[1:(n-1)])
     X <- rbind(X[1,] * (1 - rho^2)^0.5, X[2:n,] - rho * X[1:(n-1),])
     mod <- lm(y ~ X - 1)
     result <- list()
     result$coefficients <- coef(mod)
     names(result$coefficients) <- names
     summary <- summary(mod, corr = F)
     result$cov <- (summary$sigma^2) * summary$cov.unscaled
     dimnames(result$cov) <- list(names, names)
     result$sigma <- summary$sigma
     result$rho <- rho
     class(result) <- 'prais.winsten'
     result
     }

######

### Apply Prais-Winsten transformation to estimated model
PWMod <- prais.winsten.lm(Mod1)
PWMod
PWcoefs <- as.matrix(PWMod$coefficients)
xtable(PWcoefs, digits=3)
PWMod$rho
### This doesn't recover the estimates produced by the authors 

### Results from OLS & PW models
results <- cbind(OLScoefs, PWcoefs)[1:13,]
xtable(results, digits=3)



############################################
### Obtain panel-corrected standard errors

# Substitute NAs w/ zeros in the dataset, 
# as recommended in the documentation of the pcse pkg
# in the case of unbalanced panels
# (otherwise pcse issues an error message)
tsdata0 <- tsdata
tsdata0[is.na(tsdata0)] <- 0

### Simple OLS regression
Mod2 <- lm(i_entitlement ~ L.i_entitlement + D.netimmrate + D.laborcost_new +
  Dimmrate_Dlabcost + Dimmrate_Lsocial + Dimmrate_Lleft +
  L.i_socialcorp3 + L.leftgov + D.unemp + D.trade + D.gdp  + D.fdi +
  as.factor(state), data=tsdata0)
summary(Mod2)

Mod2.pcse <- pcse(Mod2, groupN=tsdata0$state, groupT=tsdata0$year, pairwise=TRUE)
summary(Mod2.pcse)
### Note that following the recommendation in the pcse documentation
# no longer produces errors,
# however this also changes the coefficients, which should not be the case
# and still doesnt recover the PCSEs in the article



###############################################


#########################################################################################
#           Replication of Tables and Figures Based on a Simple OLS Model               #
#########################################################################################



### Create function to produce expected values & 95% CI across simulations
exp.v.fn <- function(setX, simB){
	Bhat <- setX%*%simB
    ev <- apply(Bhat,1,mean)
    lb95ci <- apply(Bhat,1,quantile,.025)
    ub95ci <- apply(Bhat,1,quantile,.975)
    exp.v <- cbind(ev,lb95ci,ub95ci)      
    return(exp.v)
    }
    
### Draw simulated Beta coefficients from a multivariate normal based on previoulsy estimated OLS model
simB <- t(rmvnorm(1000, mean=Mod1$coef, sigma=summary(Mod1)$cov.unscaled))    
   
   
#############################
#     Figure 2 on p. 656    #
#############################
 

### Set covariate levels 

 ## Panel (1) in Figure 2
   # Change in integration held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_lo <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))
setX_lo[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_lo[,"D.laborcost_new"] <- quantile(tsdata[,"D.laborcost_new"], .10, na.rm=TRUE)
setX_lo[,"Dimmrate_Dlabcost"] <- setX_lo[,"D.netimmrate"]*setX_lo[,"D.laborcost_new"] 
setX_lo[,"Dimmrate_Lsocial"] <- setX_lo[,"D.netimmrate"]*setX_lo[,"L.i_socialcorp3"] 
setX_lo[,"Dimmrate_Lleft"] <- setX_lo[,"D.netimmrate"]*setX_lo[,"L.leftgov"] 


## Panel (2) in Figure 2
   # Change in integration held @ 90th pctile (= high level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_hi <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))
setX_hi[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_hi[,"D.laborcost_new"] <- quantile(tsdata[,"D.laborcost_new"], .90, na.rm=TRUE)
setX_hi[,"Dimmrate_Dlabcost"] <- setX_hi[,"D.netimmrate"]*setX_hi[,"D.laborcost_new"] 
setX_hi[,"Dimmrate_Lsocial"] <- setX_hi[,"D.netimmrate"]*setX_hi[,"L.i_socialcorp3"] 
setX_hi[,"Dimmrate_Lleft"] <- setX_hi[,"D.netimmrate"]*setX_hi[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
## Estimates for Panel (1) in Figure 2
est.lo <- exp.v.fn(setX_lo, simB)
## Estimates for Panel (2) in Figure 2
est.hi <- exp.v.fn(setX_hi, simB)

### Generate plot
pdf("AJPS2011_Figure2.pdf", 10, 6)
par(mfrow=c(1,2))
 ## Panel (1) in Figure 2
plot(0, type="n", xlim=range(-10,10), ylim=range(28,33),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo[,"D.netimmrate"], rev(setX_lo[,"D.netimmrate"])),
        y = c(est.lo[, "ub95ci"], rev(est.lo[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo[,"D.netimmrate"], est.lo[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32, "Low Level of Change \n in Integration")
legend(-5,28.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)       
  ## Panel (2) in Figure 2
plot(0, type="n", xlim=range(-10,10), ylim=range(28,33),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi[,"D.netimmrate"], rev(setX_hi[,"D.netimmrate"])),
        y = c(est.hi[,"ub95ci"], rev(est.hi[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi[,"D.netimmrate"], est.hi[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32, "High Level of Change \n in Integration")
legend(-5,28.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)
  ## Add plot title
par(mfrow=c(1,1))
mtext("FIGURE 2: Levels of Change in Integration Across the Range of Change in Immigration \n All other variables are held at their mean values \n")
dev.off()

#########  End of Figure 2 ########

 
#############################
#     Figure 3 on p. 656    #
#############################
 

### Set covariate levels 
## Panel (1) in Figure 3
   # Change in integration varies from min to max 
   # Change in immigration held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # all other covariates held @ their mean
setX_lo3 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_lo3[,"D.laborcost_new"]  <- seq(min(tsdata[,"D.laborcost_new"], na.rm=TRUE), max(tsdata[,"D.laborcost_new"], na.rm=TRUE),  len=100)
setX_lo3[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .10, na.rm=TRUE)


setX_lo3[,"Dimmrate_Dlabcost"] <- setX_lo3[,"D.netimmrate"]*setX_lo3[,"D.laborcost_new"] 
setX_lo3[,"Dimmrate_Lsocial"] <- setX_lo3[,"D.netimmrate"]*setX_lo3[,"L.i_socialcorp3"] 
setX_lo3[,"Dimmrate_Lleft"] <- setX_lo3[,"D.netimmrate"]*setX_lo3[,"L.leftgov"] 


## Panel (2) in Figure 3
   # Change in integration held @ 90th pctile (= high level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean

setX_hi3 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_hi3[,"D.laborcost_new"]  <- seq(min(tsdata[,"D.laborcost_new"], na.rm=TRUE), max(tsdata[,"D.laborcost_new"], na.rm=TRUE),  len=100)
setX_hi3[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .90, na.rm=TRUE)

setX_hi3[,"Dimmrate_Dlabcost"] <- setX_hi3[,"D.netimmrate"]*setX_hi3[,"D.laborcost_new"] 
setX_hi3[,"Dimmrate_Lsocial"] <- setX_hi3[,"D.netimmrate"]*setX_hi3[,"L.i_socialcorp3"] 
setX_hi3[,"Dimmrate_Lleft"] <- setX_hi3[,"D.netimmrate"]*setX_hi3[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
## Estimates for Panel (1) in Figure 2
est.lo3 <- exp.v.fn(setX_lo3, simB)
## Estimates for Panel (2) in Figure 2
est.hi3 <- exp.v.fn(setX_hi3, simB)

### Generate plot
pdf("AJPS2011_Figure3.pdf", 10, 6)
par(mfrow=c(1,2))

 ## Panel (1) in Figure 3
plot(0, type="n", xlim=range(-.15,.2), ylim=range(29,32.5),
     xlab ="Change in Integration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo3[,"D.laborcost_new"], rev(setX_lo3[,"D.laborcost_new"])),
        y = c(est.lo3[,"ub95ci"], rev(est.lo3[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo3[,"D.laborcost_new"], est.lo3[,"ev"], col="black", type="l", lwd=2, lty=2)
text(0, 32, "Low Level of Change \n in Immigration")
legend(-.10,29.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)
       
  ## Panel (2) in Figure 3
plot(0, type="n", xlim=range(-.15,.2), ylim=range(29,32.5),
     xlab ="Change in Integration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi3[,"D.laborcost_new"], rev(setX_hi3[,"D.laborcost_new"])),
        y = c(est.hi3[,"ub95ci"], rev(est.hi3[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi3[,"D.laborcost_new"], est.hi3[,"ev"], col="black", type="l", lwd=2, lty=2)
text(0, 32, "High Level of Change \n in Immigration")
legend(-.10,29.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)
## Add plot title
par(mfrow=c(1,1))
mtext("FIGURE 3: Levels of Change in Immigration Across the Range of Change in Integration \n All other variables are held at their mean values \n")
dev.off()
#### 
#########  End of Figure 3 ########

 
#############################
#     Figure 4 on p. 657    #
#############################

### Set covariate levels 

####  ## Panel (1) in Figure 4
   # Change in Left Gov held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_lo4 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_lo4[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_lo4[,"L.leftgov"] <- quantile(tsdata[,"L.leftgov"], .10, na.rm=TRUE)

setX_lo4[,"Dimmrate_Dlabcost"] <- setX_lo4[,"D.netimmrate"]*setX_lo4[,"D.laborcost_new"] 
setX_lo4[,"Dimmrate_Lsocial"] <- setX_lo4[,"D.netimmrate"]*setX_lo4[,"L.i_socialcorp3"] 
setX_lo4[,"Dimmrate_Lleft"] <- setX_lo4[,"D.netimmrate"]*setX_lo4[,"L.leftgov"] 


## Panel (2) in Figure 4
   # Change in Left Gov held @ 90th pctile (= high level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_hi4 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_hi4[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_hi4[,"L.leftgov"] <- quantile(tsdata[,"L.leftgov"], .90, na.rm=TRUE)

setX_hi4[,"Dimmrate_Dlabcost"] <- setX_hi4[,"D.netimmrate"]*setX_hi4[,"D.laborcost_new"] 
setX_hi4[,"Dimmrate_Lsocial"] <- setX_hi4[,"D.netimmrate"]*setX_hi4[,"L.i_socialcorp3"] 
setX_hi4[,"Dimmrate_Lleft"] <- setX_hi4[,"D.netimmrate"]*setX_hi4[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
## Estimates for Panel (1) in Figure 2
est.lo4 <- exp.v.fn(setX_lo4, simB)
## Estimates for Panel (2) in Figure 2
est.hi4 <- exp.v.fn(setX_hi4, simB)

### Generate plot
pdf("AJPS2011_Figure4.pdf", 10, 6)
par(mfrow=c(1,2))

## Panel (1) in Figure 4
plot(0, type="n", xlim=range(-10,10), ylim=range(28,33),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo4[,"D.netimmrate"], rev(setX_lo4[,"D.netimmrate"])),
        y = c(est.lo4[,"ub95ci"], rev(est.lo4[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo4[,"D.netimmrate"], est.lo4[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32.5, "Low Level of Left Party Seats")
legend(-5,28.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)
       
  ## Panel (2) in Figure 4
plot(0, type="n", xlim=range(-10,10), ylim=range(28,33),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi4[,"D.netimmrate"], rev(setX_hi4[,"D.netimmrate"])),
        y = c(est.hi4[,"ub95ci"], rev(est.hi4[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi4[,"D.netimmrate"], est.hi4[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32.5, "High Level of  Left party Seats")
legend(-5,28.5, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)

## Add plot title
par(mfrow=c(1,1))
mtext("FIGURE 4: Levels of Left Party Seats Across the Range of Change in Immigration \n All other variables are held at their mean values \n")
dev.off()
#### 
#########  End of Figure 4 ########





#############################
#     Figure 6 on p. 657    #
#############################

### Set covariate levels 


### Panel (1) in Figure 5
   # Change in left party seats varies from minimum to maximum
   # Change in immigration held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # all other covariates held @ their mean

setX_lo5 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_lo5[,"L.leftgov"]  <- seq(min(tsdata[,"L.leftgov"], na.rm=TRUE), max(tsdata[,"L.leftgov"], na.rm=TRUE),  len=100)
setX_lo5[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .10, na.rm=TRUE)

setX_lo5[,"Dimmrate_Dlabcost"] <- setX_lo5[,"D.netimmrate"]*setX_lo5[,"D.laborcost_new"] 
setX_lo5[,"Dimmrate_Lsocial"] <- setX_lo5[,"D.netimmrate"]*setX_lo5[,"L.i_socialcorp3"] 
setX_lo5[,"Dimmrate_Lleft"] <- setX_lo5[,"D.netimmrate"]*setX_lo5[,"L.leftgov"] 

### Panel (2) in Figure 5
   # Change in left party seat varies from minimum to maz
   # change in immigration held @ 90th pctile (= high level as defined by the authors on p. 655),   
   # all other covariates held @ their mean

setX_hi5 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_hi5[,"L.leftgov"]  <- seq(min(tsdata[,"L.leftgov"], na.rm=TRUE), max(tsdata[,"L.leftgov"], na.rm=TRUE),  len=100)
setX_hi5[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .90, na.rm=TRUE)

setX_hi5[,"Dimmrate_Dlabcost"] <- setX_hi5[,"D.netimmrate"]*setX_hi5[,"D.laborcost_new"] 
setX_hi5[,"Dimmrate_Lsocial"] <- setX_hi5[,"D.netimmrate"]*setX_hi5[,"L.i_socialcorp3"] 
setX_hi5[,"Dimmrate_Lleft"] <- setX_hi5[,"D.netimmrate"]*setX_hi5[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
	## Estimates for Panel (1) in Figure 2
est.lo5 <- exp.v.fn(setX_lo5, simB)

	## Estimates for Panel (2) in Figure 2
est.hi5 <- exp.v.fn(setX_hi5, simB)

### Generate plot
pdf("AJPS2011_Figure5.pdf", 10, 6)
par(mfrow=c(1,2))


#### Panel (1) in Figure 5
plot(0, type="n", xlim=range(0, 70), ylim=range(26,34),
     xlab ="Left Party Seats", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo5[,"L.leftgov"], rev(setX_lo5[,"L.leftgov"])),
        y = c(est.lo5[,"ub95ci"], rev(est.lo5[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo5[,"L.leftgov"], est.lo5[,"ev"], col="black", type="l", lwd=2, lty=2)
text(35, 32, "Low Level of Change \n in Immigration")
legend(20,28, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)
       
#### Panel (2) in Figure 5
plot(0, type="n", xlim=range(0,70), ylim=range(26,34),
     xlab ="Left Party Seats", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi5[,"L.leftgov"], rev(setX_hi5[,"L.leftgov"])),
        y = c(est.hi5[,"ub95ci"], rev(est.hi5[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi5[,"L.leftgov"], est.hi5[,"ev"], col="black", type="l", lwd=2, lty=2)
text(35, 32, "High Level of Change \n in Immigration")
legend(20,28, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)

dev.off()

#### 

#########  End of Figure 5 ########


#############################
#     Figure 6 on p. 658    #
#############################

### Set covariate levels 

####  ## Panel (1) in Figure 6
   # Change in Union Density held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_lo6 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_lo6[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_lo6[,"L.i_socialcorp3"] <- quantile(tsdata[,"L.i_socialcorp3"], .10, na.rm=TRUE)

setX_lo6[,"Dimmrate_Dlabcost"] <- setX_lo6[,"D.netimmrate"]*setX_lo6[,"D.laborcost_new"] 
setX_lo6[,"Dimmrate_Lsocial"] <- setX_lo6[,"D.netimmrate"]*setX_lo6[,"L.i_socialcorp3"] 
setX_lo6[,"Dimmrate_Lleft"] <- setX_lo6[,"D.netimmrate"]*setX_lo6[,"L.leftgov"] 


## Panel (2) in Figure 6
   # Change in Left Gov held @ 90th pctile (= high level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean
setX_hi6 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_hi6[,"D.netimmrate"]  <- seq(min(tsdata[,"D.netimmrate"], na.rm=TRUE), max(tsdata[,"D.netimmrate"], na.rm=TRUE),  len=100)
setX_hi6[,"L.i_socialcorp3"] <- quantile(tsdata[,"L.i_socialcorp3"], .90, na.rm=TRUE)

setX_hi6[,"Dimmrate_Dlabcost"] <- setX_hi6[,"D.netimmrate"]*setX_hi6[,"D.laborcost_new"] 
setX_hi6[,"Dimmrate_Lsocial"] <- setX_hi6[,"D.netimmrate"]*setX_hi6[,"L.i_socialcorp3"] 
setX_hi6[,"Dimmrate_Lleft"] <- setX_hi6[,"D.netimmrate"]*setX_hi6[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
## Estimates for Panel (1) in Figure 2
est.lo6 <- exp.v.fn(setX_lo6, simB)
## Estimates for Panel (2) in Figure 2
est.hi6 <- exp.v.fn(setX_hi6, simB)

### Generate plot
pdf("AJPS2011_Figure6.pdf", 10, 6)
par(mfrow=c(1,2))

## Panel (1) in Figure 6
plot(0, type="n", xlim=range(-10,10), ylim=range(27,34),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo6[,"D.netimmrate"], rev(setX_lo6[,"D.netimmrate"])),
        y = c(est.lo6[,"ub95ci"], rev(est.lo6[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo6[,"D.netimmrate"], est.lo6[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32.5, "Low Level of Union Density")
legend(-5,27.75, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)
       
  ## Panel (2) in Figure 6
plot(0, type="n", xlim=range(-10,10), ylim=range(27,34),
     xlab ="Change in Immmigration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi6[,"D.netimmrate"], rev(setX_hi6[,"D.netimmrate"])),
        y = c(est.hi6[,"ub95ci"], rev(est.hi6[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi6[,"D.netimmrate"], est.hi6[,"ev"], col="black", type="l", lwd=2, lty=2)
text(-1, 32.5, "High Level of Union Density")
legend(-5,27.75, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)

## Add plot title
par(mfrow=c(1,1))
mtext("FIGURE 6: Levels of Union Density Across the Range of Change in Immigration \n All other variables are held at their mean values \n")
dev.off()
#### 
#########  End of Figure 6 ########

#### 


#############################
#     Figure 7 on p. 658    #
#############################

### Set covariate levels 


 ## Panel (1) in Figure 7
   # Change in integration varies from min to max 
   # Change in Union Density held @ 10th pctile (= low level as defined by the authors on p. 655),  
   # all other covariates held @ their mean
setX_lo7 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_lo7[,"L.i_socialcorp3"]  <- seq(min(tsdata[,"L.i_socialcorp3"], na.rm=TRUE), max(tsdata[,"L.i_socialcorp3"], na.rm=TRUE),  len=100)
setX_lo7[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .10, na.rm=TRUE)


setX_lo7[,"Dimmrate_Dlabcost"] <- setX_lo7[,"D.netimmrate"]*setX_lo7[,"D.laborcost_new"] 
setX_lo7[,"Dimmrate_Lsocial"] <- setX_lo7[,"D.netimmrate"]*setX_lo7[,"L.i_socialcorp3"] 
setX_lo7[,"Dimmrate_Lleft"] <- setX_lo7[,"D.netimmrate"]*setX_lo7[,"L.leftgov"] 


## Panel (2) in Figure 2
   # Change in integration held @ 90th pctile (= high level as defined by the authors on p. 655),  
   # change in immigration varies from min to max, 
   # all other covariates held @ their mean

setX_hi7 <- t(replicate(100, apply(model.matrix(Mod1), 2, mean, na.rm=TRUE), simplify=TRUE))

setX_hi7[,"L.i_socialcorp3"]  <- seq(min(tsdata[,"L.i_socialcorp3"], na.rm=TRUE), max(tsdata[,"L.i_socialcorp3"], na.rm=TRUE),  len=100)
setX_hi7[,"D.netimmrate"] <- quantile(tsdata[,"D.netimmrate"], .90, na.rm=TRUE)

setX_hi7[,"Dimmrate_Dlabcost"] <- setX_hi7[,"D.netimmrate"]*setX_hi7[,"D.laborcost_new"] 
setX_hi7[,"Dimmrate_Lsocial"] <- setX_hi7[,"D.netimmrate"]*setX_hi7[,"L.i_socialcorp3"] 
setX_hi7[,"Dimmrate_Lleft"] <- setX_hi7[,"D.netimmrate"]*setX_hi7[,"L.leftgov"] 


### Compute predicted values based on simulated coefficients and specified covariate levels
## Estimates for Panel (1) in Figure 2
est.lo7 <- exp.v.fn(setX_lo7, simB)
## Estimates for Panel (2) in Figure 2
est.hi7 <- exp.v.fn(setX_hi7, simB)

### Generate plot
pdf("AJPS2011_Figure7.pdf", 10, 6)
par(mfrow=c(1,2))


 ## Panel (1) in Figure 3
plot(0, type="n", xlim=range(.1,.9), ylim=range(27,33),
     xlab ="Change in Integration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_lo7[,"L.i_socialcorp3"], rev(setX_lo7[,"L.i_socialcorp3"])),
        y = c(est.lo7[,"ub95ci"], rev(est.lo7[,"lb95ci"])), col="grey55", border=NA)
lines(setX_lo7[,"L.i_socialcorp3"], est.lo7[,"ev"], col="black", type="l", lwd=2, lty=2)
text(.45, 32, "Low Level of Union Density")
legend(.15,28, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey55", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey55", NA), pt.cex = 3)
       
  ## Panel (2) in Figure 2
plot(0, type="n", xlim=range(.1,.9), ylim=range(27,33),
     xlab ="Change in Integration", ylab="Predicted Level of Unemployment Entitlement")
polygon(x = c(setX_hi7[,"L.i_socialcorp3"], rev(setX_hi7[,"L.i_socialcorp3"])),
        y = c(est.hi7[,"ub95ci"], rev(est.hi7[,"lb95ci"])), col="grey90", border=NA)
lines(setX_hi7[,"L.i_socialcorp3"], est.hi7[,"ev"], col="black", type="l", lwd=2, lty=2)
text(0.45, 32, "High Level of Union Density")

legend(.15,28, legend = c("95% Confidence Interval", "Prediction"), cex=.8,
       col = c("grey90", "black"), lty = c(0, 2), lwd = c(0, 2),
       pch = c(22, NA), pt.bg = c("grey90", NA), pt.cex = 3)

dev.off()
#### 

